Load network and scent data
From Tables S1 and S2 in: Kantsa, A., R.A. Raguso, A.G. Dyer, J.M. Olesen, T. Tscheulin, and T. Petanidou. 2018. Disentangling the role of floral sensory stimuli in pollination networks. Nature Communications 9: 1041. doi:[10.1038/s41467-018-03448-w](https://doi.org/10.1038/s41467-018-03448-w)
#load network
net <- read.delim("Kantsa_etal_Supp_Data1.csv", skip = 2)
rownames(net) <- net[,1]
net[,1] <- NULL
net <- as.matrix(t(net))
polls <- colnames(net)
snet <- sortweb(net) #sort by row & column totals
#load floral scents
scent.file <- "Kantsa_2018_Supp_Data2.xlsx"
plants <- getSheetNames(scent.file)
scentlist <- lapply(2:length(plants), read.xlsx, xlsxFile=scent.file)
plants <- plants[-1]
scent <- rbindlist(scentlist, fill=T, idcol="plant")[,-4]
colnames(scent)[2:3] <- c("cmpd","mean")
scent$plant <- factor(plants[scent$plant])
scent$cmpd <- factor(tolower(scent$cmpd))
cmpds <- levels(scent$cmpd)
scent.net <- dcast(scent, plant ~ cmpd, value.var="mean", fun.aggregate=mean, fill=0) #long to wide
scent.net <- as.data.frame(scent.net)
rownames(scent.net) <- scent.net[,1]
scent.net[,c(1,380)] <- NULL #weird NA column has one entry
scent.net <- as.matrix(scent.net)
scent.net.short <- scent.net[,colSums(scent.net)>0.1]Plant phylogeny
#plants.class <- classification(plants, db="gbif")
#save(plants.class, file="plantsclass.Rdata")
load("plantsclass.Rdata")
plants.classtree <- class2tree(plants.class) #taxonomic tree
plants.tree <-phylomatic(taxa=names(plants.class), get = 'POST', storedtree="zanne2014")
plants.tree$tip.label[order(plants.tree$tip.label)] <- plants
par(mar=c(0,0,0,0))
cp <- cophylo(plants.tree, plants.classtree$phylo) Rotating nodes to optimize matching...
Done.plot(cp)
Pollinator phylogeny
polls <- gsub("Bug", "Hemiptera", polls, fixed=T)
polls <- gsub("Auchenorrhyncha", "Hemiptera", polls, fixed=T)
polls <- gsub("Fly", "Diptera", polls, fixed=T)
polls <- gsub("Beetle", "Coleoptera", polls, fixed=T)
#polls.class <- classification(polls, db="gbif")
#save(polls.class, file="pollsclass.Rdata")
load("pollsclass.Rdata")
polls.gbif <- factor(sapply(polls.class, function(x) tail(x$name, 1)))
net.lump <- aggregate(.~ polls.gbif, as.data.frame(t(net)), sum)
rownames(net.lump) <- net.lump[,1]
net.lump[,1] <- NULL
net.lump <- as.matrix(t(net.lump))
rownames(net.lump) <- sort(plants.tree$tip.label)
polls.classtree <- class2tree(unique(polls.class)) #taxonomic tree
polls.matched <- tnrs_match_names(lapply(polls.class, function(x) tail(x$name, 1)), context_name = "Animals") Warning: Some names were duplicated: 'andrena', 'asilidae', 'coleoptera',
'hemiptera', 'buprestidae', 'buprestidae', 'calliphoridae', 'cerambycidae',
'chrysididae', 'crabronidae', 'eumenidae', 'eumenidae', 'eumenidae',
'diptera', 'diptera', 'diptera', 'diptera', 'diptera', 'hoplitis',
'hoplitis', 'lasioglossum', 'noctuidae', 'oedemera', 'osmia', 'osmia'. Warning in check_tnrs(res): chrysura circe are not matchedpolls.ids <- ott_id(polls.matched)
polls.ids <- polls.ids[!polls.ids %in% c(968359, 707236, 707875, 63881, 458953, 42699, 742352, 340559, 672497, 34898, 332969,356987)]
polls.tree <- tol_induced_subtree(ott_ids=polls.ids)
Progress [-----------------------------------] 0/7 ( 0%) ?s
Progress [==================================] 7/7 (100%) 0s
Warning in collapse_singles(tr): Dropping singleton nodes with labels:
Colletes ott115493, Crabronidae ott372234, Scolia ott881314, Coleoptera
ott865243, Diptera ott661378, Paragus ott973598, Hemiptera ott603650par(mar=c(0,0,0,0))
plot(polls.classtree, cex=0.5)
plot(polls.tree, cex=0.5)
Plant-pollinator network
Cophylogeny
mnet.lump <- melt(net.lump)
mnet.lump <- mnet.lump[mnet.lump$value>0,]
cophyloplot(plants.tree, polls.classtree$phylo, mnet.lump[,1:2], lwd=sqrt(mnet.lump[,3])/4, show.tip.label=F, col='steelblue', length.line=0, gap=-20, space=60)
Pollinator groups
polls.cut <- dendextend::cutree(polls.classtree$phylo,5)
polls.order <- setNames(factor(sapply(unique(polls.class), function(x) x$name[4])), names(polls.class[!duplicated(polls.class)]))
polls.order.all <- setNames(factor(sapply(polls.class, function(x) x$name[4])),names(polls.class))
#table(polls.cut, polls.order)
net.grp <- aggregate(.~ polls.order.all, as.data.frame(t(net)), sum)
rownames(net.grp) <- net.grp[,1]
net.grp[,1] <- NULL
net.grp <- as.matrix(t(net.grp))
rownames(net.grp) <- sort(plants.tree$tip.label)
plants.majorpoll <- setNames(factor(colnames(net.grp)[apply(net.grp, 1, which.max)]), rownames(net.grp))
pcols <- setNames(brewer.pal(4,"Set1"), levels(plants.majorpoll))
plants.polltrees<-make.simmap(plants.tree,plants.majorpoll,model="ARD",nsim=100) make.simmap is sampling character histories conditioned on the transition matrix
Q =
Coleoptera Diptera Hymenoptera Lepidoptera
Coleoptera -0.008667015 0.0006377689 0.008029247 0.00000000
Diptera 0.017012461 -0.0224322446 0.005419784 0.00000000
Hymenoptera 0.000000000 0.0000000000 -0.071181044 0.07118104
Lepidoptera 0.078649350 0.0796271485 1.270049241 -1.42832574
(estimated using likelihood);
and (mean) root node prior probabilities
pi =
Coleoptera Diptera Hymenoptera Lepidoptera
0.25 0.25 0.25 0.25 Done.plants.polltree<-summary(plants.polltrees,plot=FALSE)
plot(plants.polltree,colors=pcols,fsize=0.8,cex=c(0,0.3))
add.simmap.legend(colors=pcols,x=0.9*par()$usr[1],
y=0.2*par()$usr[4],prompt=FALSE,fsize=0.9)
pcols.all <- setNames(c(pcols[1:2],brewer.pal(5,"Set1")[5],pcols[3:4]), levels(polls.order.all))Heatmap
net.cut <- net[rowSums(net)>0,]
heatmaply(net.cut, scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
distfun=vegdist,
hclustfun=function(x) hclust(x, method="mcquitty"),
row_side_colors = data.frame(P=plants.majorpoll[rowSums(net)>0]),
row_side_palette=function(n) pcols,
col_side_colors=data.frame(PollinatorOrder=polls.order.all),
col_side_palette=function(n) pcols.all) Warning in heatmaply.heatmapr(hm, colors = colors, limits = limits,
scale_fill_gradient_fun = scale_fill_gradient_fun, : The hover text for
col_side_colors is currently not implemented (due to an issue in plotly).
We hope this would get resolved in future releases.Matrix plot
#interaction matrix plot
visweb(net, type="nested", circles=T, boxes=F, circle.max=1.2)
Web plot
#interaction web plot
plotweb(sqrt(net), empty=F, text.rot=90, method="cca", col.high=pcols.all[polls.order.all], col.low=pcols[plants.majorpoll], col.interaction=pcols.all[polls.order.all], bor.col.interaction=pcols.all[polls.order.all], bor.col.high=pcols.all[polls.order.all], bor.col.low=pcols[plants.majorpoll])
Network topography
networklevel(net) connectance web asymmetry
0.06312657 0.63106796
links per species number of compartments
1.95631068 2.00000000
compartment diversity cluster coefficient
1.07901412 0.02631579
nestedness weighted nestedness
4.47259022 0.58578362
weighted NODF interaction strength asymmetry
10.80348739 0.20372363
specialisation asymmetry linkage density
-0.04424002 3.57255189
weighted connectance Fisher alpha
0.01734248 89.61252715
Shannon diversity interaction evenness
3.27629195 0.37393976
Alatalo interaction evenness H2
0.25285628 0.57919894
number.of.species.HL number.of.species.LL
168.00000000 38.00000000
mean.number.of.shared.partners.HL mean.number.of.shared.partners.LL
0.34944397 1.18776671
cluster.coefficient.HL cluster.coefficient.LL
0.34062967 0.21733898
weighted.cluster.coefficient.HL weighted.cluster.coefficient.LL
0.72421977 0.31277076
niche.overlap.HL niche.overlap.LL
0.12366249 0.14889538
togetherness.HL togetherness.LL
0.03019814 0.02302573
C.score.HL C.score.LL
0.75923950 0.74100046
V.ratio.HL V.ratio.LL
3.23242952 17.02718828
discrepancy.HL discrepancy.LL
248.00000000 250.00000000
extinction.slope.HL extinction.slope.LL
4.79325165 1.62192500
robustness.HL robustness.LL
0.80798877 0.61786266
functional.complementarity.HL functional.complementarity.LL
8044.31456516 7651.23911634
partner.diversity.HL partner.diversity.LL
0.94103312 1.28483306
generality.HL vulnerability.LL
2.67434342 4.47076037Scent network
Heatmap
heatmaply(sqrt(scent.net), scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
distfun=vegdist,
hclustfun=function(x) hclust(x, method="mcquitty"),
row_side_colors = data.frame(P=plants.majorpoll),
row_side_palette=function(n) pcols)Web plot
#interaction web plot
plotweb(sqrt(scent.net), text.rot=90, method="cca", col.high=pal[2], bor.col.high=pal[2], col.low=pcols[plants.majorpoll], bor.col.low=pcols[plants.majorpoll], col.interaction=pal[5], bor.col.interaction=pal[5])
Plant-scent-pollinator network
NMDS
scent.mds <- metaMDS(sqrt(scent.net), autotransform = F, trace=F)
#scent.op <- ordiplot(scent.mds, type="n")
#points(scent.op, what="species", pch=3, col="grey50")
#points(scent.op, what="sites", cex=1.5, pch=19, col=pcols[plants.majorpoll])
#text(scent.op, what="sites", cex=0.8, col=pcols[plants.majorpoll])Phylochemospace
scent.mds.taxa <- scent.mds$points
rownames(scent.mds.taxa) <- sort(plants.tree$tip.label)
rownames(scent.mds.taxa) <- gsub(" ", "_", rownames(scent.mds.taxa), fixed=T)
plants.tree$tip.label <- gsub(" ", "_", plants.tree$tip.label, fixed=T)
#plants.tree$phylo$edge.length <- plants.tree$phylo$edge.length+0.00001
#Do generalists cluster?
plants.numlinks <- rowSums(decostand(net.lump, "pa"))
tip.pcols <- setNames(colorRampPalette(c("white", "darkblue"))(67)[plants.numlinks+1],
sort(plants.tree$tip.label))
node.pcols<- setNames(c(tip.pcols[plants.tree$tip.label], rep("black",plants.tree$Nnode)),
1:(length(plants.tree$tip)+plants.tree$Nnode))
phylomorphospace(plants.tree, scent.mds.taxa, control=list(col.node=node.pcols), node.size=c(0,2), label="off", lwd=1)
text(scent.mds.taxa[,1], scent.mds.taxa[,2], rownames(scent.mds.taxa), cex=0.8, offset=0.5, xpd=T, col=tip.pcols)
#Do plant major pollinators cluster?
tip.pcols <- setNames(pcols[plants.majorpoll],
sort(plants.tree$tip.label))
node.pcols<- setNames(c(tip.pcols[plants.tree$tip.label], rep("black",plants.tree$Nnode)),
1:(length(plants.tree$tip)+plants.tree$Nnode))
phylomorphospace(plants.tree, scent.mds.taxa, control=list(col.node=node.pcols), node.size=c(0,2), label="off", lwd=1)
text(scent.mds.taxa[,1], scent.mds.taxa[,2], rownames(scent.mds.taxa), cex=0.8, offset=0.5, xpd=T, col=tip.pcols)
library(vegan3d)
orditree3d(scent.mds.taxa, as.hclust(force.ultrametric(multi2di(plants.tree))), lwd=2, pch=19, col=tip.pcols)
phenogram(plants.tree, setNames(plants.numlinks,sort(plants.tree$tip.label))) Optimizing the positions of the tip labels...
png(file="ppm-%04d.png",width=600,height=600,res=120)
par(mar=c(2.1,2.1,1.1,1.1))
project.phylomorphospace(tree=plants.tree, X=scent.mds.taxa, xlab="", ylab="", node.size=c(0,0), lwd=1, direction="to", nsteps=100, fsize=0.6)
dev.off()
system("convert -delay 10 -loop 2 *.png phylochemo.gif")
file.remove(list.files(pattern=".png"))Phylogeny scent bar chart
ccol = sample(rainbow(ncol(scent.net)))
par(mfrow=c(1,2))
par(mar=c(0.8,0,0.8,0))
plot(plants.tree, cex=1, tip.color = tip.pcols)
par(mar=c(0,0,0,0))
barplot(t(scent.net), col=ccol, names.arg=rep("",nrow(scent.net)), horiz=TRUE) 
par(mfrow=c(1,2))
par(mar=c(0.8,0,0.8,0))
plot(plants.tree, cex=1, tip.color = tip.pcols)
par(mar=c(0,0,0,0))
barplot(t(decostand(scent.net, method="tot")), col=ccol, names.arg=rep("",nrow(scent.net)), horiz=TRUE) 
Indicator compounds
library(indicspecies)
mp <- multipatt(as.data.frame(scent.net), plants.majorpoll, control=how(nperm=999))
summary(mp)
Multilevel pattern analysis
---------------------------
Association function: IndVal.g
Significance level (alpha): 0.05
Total number of species: 378
Selected number of species: 2
Number of species associated to 1 group: 2
Number of species associated to 2 groups: 0
Number of species associated to 3 groups: 0
List of species associated to each combination:
Group Lepidoptera #sps. 2
stat p.value
un 68, 67, 81, 54, 41, 95 1.000 0.027 *
3-carene 0.986 0.014 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(paco)
# start the paco procedure
poll_phy <- polls.classtree$phylo
pla_phy <- as.phylo(hclust(vegdist(scent.net), "mcquitty"))
int <- net.lump
H <- cophenetic(polls.classtree$phylo)
P <- cophenetic(as.phylo(hclust(vegdist(scent.net), "mcquitty")))
D <- prepare_paco_data(H, P, net.lump) Warning in prepare_paco_data(H, P, net.lump): The HP matrix should have
hosts in rows. It has been translated.D <- add_pcoord(D, correction='cailliez')
# now we are ready for cophylogenetic analysis
#D <- PACo(D, nperm=1000, seed=13, method='quasiswap', symmetric=TRUE)
#save(D, file="D.Rdata")
load("D.Rdata")
# and to investigate the contribution of individual links
res <- residuals_paco(D$proc)
# to visualise this we use the ape function cophyloplot weighted by interaction contribution
# first we must make a list out of the interaction matrix
assoc <- data.frame(pol=rownames(net.lump)[which(net.lump>0, arr.ind=TRUE)[,'row']], pla=colnames(net.lump)[which(net.lump>0, arr.ind=TRUE)[,'col']])
# to weight the interactions we use the cophylogenetic contribution transformed to best show
# the differences graphically
weight <- (res^-2)/50
mnet.lump <- melt(net.lump) #replaces assoc
mnet.lump <- mnet.lump[mnet.lump$value>0,]
#cophylo with interaction strength
cophyloplot(pla_phy, poll_phy, mnet.lump[,1:2], lwd=sqrt(mnet.lump[,3])/4, show.tip.label=F, col='steelblue', length.line=0, gap=-20, space=60)
#cophylo with residuals
cophyloplot(pla_phy, poll_phy, mnet.lump[,1:2], lwd=weight/8, show.tip.label=F, col='steelblue', length.line=0, gap=-20, space=60)
Do generalist plants have higher diversity of scent compounds?
plants.numlinks <- rowSums(decostand(net.lump, "pa"))
polls.classtree <- compute.brlen(polls.classtree$phylo,1)
library(picante) Loading required package: nlme
Attaching package: 'nlme' The following object is masked from 'package:sna':
gapplyplants.pdlinks <- pd(net.lump, polls.classtree)$PD
plants.numcompounds <- rowSums(decostand(scent.net, "pa"))
plants.sumcompounds <- rowSums(scent.net)
plants.shannoncompounds <- diversity(scent.net, index="shannon")
plants.shannonlinks <- diversity(net.lump, index="shannon")
plot(plants.numcompounds~plants.numlinks)
plot(plants.numcompounds~plants.pdlinks)
plot(plants.shannoncompounds~plants.pdlinks)
plot(plants.shannoncompounds~plants.shannonlinks)
plot(sqrt(plants.sumcompounds)~plants.pdlinks)
plot(plants.shannoncompounds~plants.majorpoll)
plot(plants.numcompounds~plants.majorpoll)
plot(sqrt(plants.sumcompounds)~plants.majorpoll)
qplot(plants.pdlinks, plants.numcompounds) + geom_smooth(method="lm")
summary(lm(plants.numcompounds[-9]~plants.pdlinks[-9]))
Call:
lm(formula = plants.numcompounds[-9] ~ plants.pdlinks[-9])
Residuals:
Min 1Q Median 3Q Max
-19.665 -10.691 -1.845 9.243 30.240
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 20.76015 4.15181 5.000 1.33e-05 ***
plants.pdlinks[-9] 0.02646 0.07241 0.365 0.717
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 13.82 on 38 degrees of freedom
Multiple R-squared: 0.003502, Adjusted R-squared: -0.02272
F-statistic: 0.1335 on 1 and 38 DF, p-value: 0.7168data(phylocom)
pd(phylocom$sample, phylocom$phylo) PD SR
clump1 16 8
clump2a 17 8
clump2b 18 8
clump4 22 8
even 30 8
random 27 8